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Abstract 

We discuss the problem of two particles interacting via short-range interactions within a 
harmonic-oscillator trap. The interactions are organized according to their number of deriva- 
tives and defined in truncated model spaces made from a bound-state basis. Leading-order (LO) 
interactions are iterated to all orders, while corrections are treated in perturbation theory. We 
show explicitly that next-to-LO and next-to-next-to-LO interactions improve convergence as the 
model space increases. In the large-model-space limit we regain results from a pseudopotential. 
Arbitrary scattering lengths are considered, as well as a generalization to include the non- vanishing 
range of the interaction. 
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I. INTRODUCTION 



Much experimental progress has recently been achieved in the area of ultra-cold atoms. 
It is now possible [H to make systems of atoms confined in optical lattices formed by laser 
beams and at the same time control the strength of the two-body interaction with magnetic 
fields. The atoms are cooled down to extremely low temperatures and, in the limit of low 
tunneling, each lattice site may be regarded as a harmonic-oscillator (HO) well, which is 
independent from the others and contains only a few atoms. Using a Feshbach resonance, 
the two-body interaction can be fine-tuned so that the S'-wave scattering length 02 is made 
much larger than the range R of the interaction. By doing so, the theoretical problem of a 
system of trapped particles interacting via short-range forces can be realized experimentally. 

Such atomic systems present remarkable similarities to some nuclear physics systems. 
Indeed, the two-nucleon scattering length is much larger than the range of the nuclear force, 
set by the pion mass. In such a situation, the physical properties of few-body systems are 
universal, that is, they depend mainly on 02, not on the details of the two-body interaction. 
At sufficiently low energies, an effective field theory (EFT) has been formulated which turns 
this separation of scales into an expansion in powers of R/a2 j2|- Except for isospin, which 
does play a role in the relative relevance of few-body forces, the version of this EFT used 
in nuclear physics jsj is formally indistinguishable from the theory describing atomic sys- 
tems (Nevertheless, the underlying theories for the two cases are very different.) As a 
consequence, atomic systems characterized by large scattering lengths can be studied with 
techniques developed in nuclear physics and, conversely, provide an excellent testing ground 
for few- and many-body methods that can be further applied, with little or no change, to 
the description of nuclear systems at low energies. 

The no-core shell model (NCSM) is one of the most fiexible ah initio methods used to ob- 
tain the solution to the non-relativistic Schrodinger equation for many-nucleon systems jsf. 
Currently, it is the only such method able to reach medium-mass nuclei with no restrictions 
to closed or nearly-closed shells and at the same time to handle local and non-local interac- 
tions on the same footing. It uses a discrete single-particle basis — typically a HO basis — 
and, being a numerical method, relies on a suitable truncation of the model space accessible 
to nucleons — in the form of a maximum number Nmax of accessible shells above the mini- 
mum configuration. This requires the use of effective interactions (not only among nucleons 
but also among nucleons and external probes such as photons, when they are present) to 
account for effects left out by the truncation to a finite model space. The method of choice 
has been the construction of effective operators via unitary transformations. This involves 
the use of an approximation, the so-called "cluster approximation", which is not a priori 
controlled, and which poses challenges for the description of low- momentum observables jsj. 
An alternative is to use EFT to construct effective interactions that are consistent with the 
underlying theory of QCD directly in the model spaces where many-nucleon calculations are 
carried out 0. 

In this paper, we consider the problem of two trapped atoms with large scattering length, 
1 02 1 ^ -R, from the perspective of EFT interactions solved with the NCSM. The main 
simplification with respect to the nuclear case is that the trapping potential provides a 
natural single-particle basis; its HO length h does not need to be removed at the end of the 
calculation because it represents the long-distance physics of the trap. As long a.sh^ R, the 
trapped system should still exhibit universal behavior, although for 6 < |a2| it is significantly 
different from that of the untrapped system. 
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We use EFT techniques to expand the interaction between particles as a series of contact 
interactions with an increasing number of derivatives, and calculate energy levels inside the 
trap by explicitly solving the Schrodinger equation. We extend here the work initiated in Ref. 
[sl, in which only the leading-order (LO) interaction was investigated. We include explicitly 
corrections up to the next-to- next-to-leading-order (N^LO), and show how they provide 
improved convergence as N^ax increases, at least as long as they are treated in perturbation 
theory. Finite and infinite values of the scattering length a2 are considered, and we show 
that in the Nmax oo limit our LO results converge to the levels of a pseudopotential 



with the same scattering length, obtained by Busch et al. [10[. Starting at NLO a finite 
value of the effective range r2 is allowed as well, and at higher orders other effective-range 
expansion parameters could be included similarly. The Nmax — > oo limit of our subleading- 



order EFT provides a derivation of the generalized Busch et al. relation [11H13| . With the 
two-body system thus understood, we can use our method to calculate the energies of larger 
trapped systems fl3 |. 

Other approaches to the same problem of course exist in the literature. They are fre- 



quently based on a specific form for the interparticle potential (see, for example, Ref. 
Closest in spirit to ours is probably the approach [l5| where an effective short-range inter- 
action is fitted to several levels of the pseudopotential at unitarity, but diagonalized exactly. 
What distinguish our framework are its i) systematic character, in the form of a controlled 
expansion; and ii) generality, since no assumptions about the form of the short-range inter- 
actions are needed. These features allow the same method to apply across fields. Of course, 
the low-energy observables themselves should be in agreement among different approaches, 
as long as they are calculated properly. 

The paper is organized as follows. We first show in Sec. [IT] the general formalism of our 
approach, including a detailed description of the treatment beyond leading order. We illus- 
trate the method in different situations (finite and infinite scattering length, negligible and 
non- vanishing effective range, etc.) in Sec. Illli We conclude and discuss future applications 
in Sec. IIVI Appendices [K \ [B | and O provide some of the details omitted in the main text. 



II. GENERAL CONSIDERATIONS 

We consider a non-relativistic system of two particles of reduced mass /i that interact 
with each other in the S wave. For definiteness, we think of two-component fermions, which 
support a single S channel and interact also in P and higher waves, but the framework can 
be straightforwardly applied to bosons and other fermions. In free space, the properties of 
this system can be characterized by scattering phase shifts for each partial wave of angular 
momentum I at relative on-shell momentum k, Si{k). When the relative momentum is much 
smaller than the inverse of the range R of the interparticle interaction, k <^ 1/R (we use 
units in which ^ = 1 and c = 1), the phase shifts are given by the effective range expansion 
(ERE); for example, in the S wave, 

k cot 6o{k) = -- + ]-r2k^ + -P2k^ + . . . , (1) 
a2 2 4 

where 02, r2, P2, • • • are, respectively, the scattering length, effective range, shape parameter, 
and higher ERE parameters not shown explicitly. A similar expansion exists for P and higher 
partial waves. 
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Generically, the sizes of ERE parameters are set by i?, for example |r2| ~ -R. The ERE ([T]) 
is an expansion in powers of kR. Given a desired precision, the ERE can be truncated and 
the system described by a finite number of parameters. Potentials that generate the same 
values for this finite number of ERE parameters cannot be distinguished at this precision 
level: they all generate the same wavefunction for distances beyond the range of the force, 
r>R. Such potentials are said to belong to the same universality class. Most interesting 
is the case where the depth of the interaction potential is fine-tuned so that an 5'-wave 
bound state is near threshold and the associated scattering length is large, \a2\ ':$> R. Then 
the physics of the bound state is largely independent of the potential; for example, for 
R r <^ a2 the wavefunction of a state of energy E is 



7/;(r) oc 3 <j 1 - [1 - i2a2r2E + . . .] ^ + C ( — 




(2) 



Physics in the unitarity hmit, \a2\ oo, is controlled by the zero-energy bound state and 
exhibits a higher degree of universality. 

The ERE description is useful because it is model independent. However, it does not 
directly provide a basis for the study of many-body systems. This can be accomplished 
using the ideas of EFT jsi 0]. Since at low momenta the details of the interaction cannot 
be resolved, we can expand the interparticle interaction V as a. Taylor series in momentum 
space, as done in App. \^ in coordinate space. 



V{f', r) = Co6{f')6{r) - C2 { [V '^S{f')] 6{r) + 6{f') [V^6{r)] } 

+C4 {[V''<5(f')] 6{r) + 6{f') [V%f)] + 2 [v'^6{f')] [v^6{r)] } 



(3) 



where Co, C2, and C4 are parameters, and ". . . " denote interactions that contribute at 
higher orders. Because the interactions are singular, an ultraviolet (UV) cutoff A has to be 
introduced to solve the Schrodinger equation. In order for observables to be independent 
of A (renormalization-group invariance), the parameters Ci have to depend on A. This is 
merely a consequence of the fact that the UV cutoff is an arbitrary separation between the 
short-range dynamics included explicitly in the dynamics (through high virtual momenta) 
versus that included implicitly in the potential (through its parameters). 

Contributions to observables obtained from V ([3]) can also be organized in powers of kR. 
It has, in fact, been shown ^ that in the two-body sector the EFT expansion reproduces 
the ERE ([1]) at each power of kR. In the generic situation, one can simply treat the whole 
potential in perturbation theory. When \a2\ S> -R, however, the Cq term in Eq. (|3]) needs 
to be solved exactly, while the remaining terms can still be accounted for in perturbation 
theory j2|. These higher-order terms represent range and subtler effects in the S wave, and 
P and higher waves. 

Here we are interested in the large-scattering-length scenario, when the two particles are 
trapped in an HO potential of frequency u. The HO introduces a third scale, the length 
b = \j J^UjJ. In the relative frame, the Hamiltonian for this system reads 



E 



2 



-6^2 + 



62 



(4) 



As long a.s b ^ R, details of the interparticle potential remain irrelevant and universality is 
not destroyed. In the following we consider various values for the ratio 6/02, which are, as 
discussed in Sec. [U of interest in both atomic and nuclear physics. 
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We want to set up a perturbative approach; we thus write the Hamiltonian for the 
relative motion as 

H = + V^'^ + ^(2) + . . . , (5) 
with the energy and wavefunction decomposed accordingly, 



1^(0)) + + 1^1^)) + ... 



■,(1) 



(2)\ 



and 



E = + + + . . . 



(6) 



(7) 



The superscript '•"^ corresponds to the order in perturbation theory of the different terms. 
Here we order interactions according to the power counting of Ref. 0] and for convenience 
we split the parameters Cj in Eq. among different orders. 
The leading-order (LO) Hamiltonian is 



u 
2 



2t-72 



r 

62 



(8) 



and the corresponding wavefunction ip^'^^r) is the solution of the Schrodinger equation 

(i7{0)_^(0))^(0)(^-.) ^ (9) 

The next-to-leading-order (NLO) correction to the potential is 

= C^^5(r) - Ci^^ { [V^S{r}] + 2 [V6{r)] ■ V + 26{r}V^} , (10) 

and the first-order corrections to the energy, E^^\ and to the wavefunction, %jj^^'>{f), are 
obtained in first-order perturbation theory. That is, they are such that 

The next-to-next-to-leading-order (N^LO) correction to the potential V^'^^ is given by 
= cS'^5(f) - ^ {[V'5(r)] + 2 [W(f)] ■ V + 25(f) V'} 



+Cf ^ { [V'(5(f)] + 4 [vV^6{f 



■ V + 4 



[VW(r 



■vv 



+4 



W(f)] ■ vv' + 25(r)V^} . (12) 
The corrections -E^') and ijjC^^l^r) are obtained from 

which just means perturbation theory to first order in V^'^^ and to second order in V^^\ 
Extension to higher orders is straightforward. 

We work on a basis of HO wavefunctions (f) with energies E^i = {2n + / + 3/2)u = 
{N + 3/2)0;. Useful properties of these wavefunctions are summarized in App. [Bl In the 
HO basis the singularity of the potential Q can be tamed by imposing a maximum number 
of shells Nmax I2I , which corresponds to a UV momentum cutoff 



A= - A/2iV^,, + 3. 



(14) 



We therefore expand the wavefunction ([6]) in the HO basis in the finite space, 

max I 

^^''\r)= Y: E ci:L<PniUrl, (15) 

n,l=0 m=—l 

where Cnim are coefficients to be determined. 

Since the two-body potentials up to N^LO only support S waves, to this order the eigen- 
functions with / > are simply the HO wavefunctions (pnimir) with eigenvalues {2n+l+3/2)uj 
up to the energy {N^ax + S/2)u. For / = 0, on the other hand, the levels are affected by 
the interparticle potential. Denoting by Umax the maximum value of the radial quantum 
number (that is, rimax is the largest integer smaller than, or equal to, Nmaxf^) and omitting 
the labels / = and m = 0, the S wavefunction can be written as 

'^max 

^t\r) = E c^rVn(r) (16) 

n=0 

in terms of the HO S wavefunctions 

0„(r) = 7r-3/4r3/2 (0)]"'^%-^'/262 4V2) (ryb') , (17) 

where L^"^ is the generalized Laguerre polynomial. 

The resulting energies ([7]) will depend on Nmax as well as u, E = E{Nmax, w). Since Nmax 
is arbitrary, we want the energies not to depend sensitively on Nmax- This cannot be achieved 
in general, but it can for the shallow levels of interest — that is, those with E < 0(l/2/xi?^), 
which are dominated by physics at distances r>R. As we show in the following, this is 

accomplished by allowing the C^'^^ to depend on both N^ax and oo, C^'^^ = C\''\Nmax,^)- 
Nevertheless, at any order a residual N^ax dependence introduces an error in the calculation 
of shallow levels, which should be proportional to powers of 1/A. At the end of the calculation 
we want to take N^ax sufficiently large, A> 1/i?, so that this error is not larger than the 
error proportional to powers of R stemming from the truncation of Eq. 



A. LO renormalization 

The physics at LO is obtained by diagonalizing the two-body Hamiltonian ([H]). The 
approach is similar to the treatment in a free-particle basis [2[, with the difference that we 
work here only with bound states, which naturally are closely related to HO states because 
of the presence of the trap. Renormalization of the interaction at LO has already been 
discussed in Ref. [8|, but for completeness we repeat the derivation here. 

We start with the Schrodinger equation ([9]) for the wavefunction of the two-fermion 
system, 'ip'^^\r). One can show, by inserting Eq. f|T6l) into Eq. ([9]) and projecting on a HO 
basis state (pnif), that the expansion coefficients are given by 



0n(O) 



- (2n + 3/2)a;' ^^^^ 

where k'-'^^ = Cg^V'^HO) is itself a combination of the unknown coefficients. Direct substi- 
tution of the expansion coefficients back into Eq. (fT6!) yields 
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From the consistency condition at r = we obtain the relation 
to satisfy: 



1 



fx 



L(V2)(0) 



'^=0 2uj v"" ^ 4/ 



2/i |r(r^ 



max ~ 2 



ir'^b 1 r (rir 



1 + R\nr, 



^{0)( 



CO 



2uj 



that an energy E^^^ has 



TT 



r 



E(°)(a;) ' 
2uj 



2w 



(?o) 



where 



R{m,- 



l-2e 



8(m + 1) (m + I - f) 

^ / 3 7 

3F2 l,m+ -,m+ - 



-: m + 2, m 



11 
T 



-;1 



(21) 



in terms of the generalized hypergeometric function 3F2. (In Eq. f l20|) we used Eqs. flB16P 
and (lB17p .) Finally, the constant k,^^^ is fixed by the chosen normalization of iljll^\r). (It 



can be calculated using Eq. (IBlSp .) Note that this constant is energy dependent; when 
necessary we denote it by fi^^(o)- 

The rhs of Eq. f l2U]) clearly depends on n.max and 00, and so does Cq°^ In a given model 
space, the coupling constant C^^ can be fixed to reproduce one observable, in this case one 
energy of the two-body system. If, to be definite, we take that energy as the ground-state 
energy in the trap. 



(0) 



^o(^), 



then CQ'\n„ 



, u) is determined from 



\''"n 




R 



2uj 



TT 



Eojuj) 
2u) 



2u} 



(22) 



.(23) 



While one energy in the spectrum is fixed in all model spaces, the rest of the energy spectrum 
runs with the model space: the remaining energies i?f>\ = El'^^{nmax,co) satisfy Eq. ( 120|1 
and in general depend not only on u but also on rimax- These energies should converge 
as Umax — 00 to finite values -E'f>\(oo, w). However, once A in Eq. iHM exceeds 1/R, the 
theoretical errors are dominated by the physics of the effective range r2, which was left out 
of LO. 



B. Renormalization beyond LO 

In the two-nucleon system, the effective range is much smaller than the scattering length, 
but is finite. In an atomic system near a Feshbach resonance, the range is usually neglected, 
although it should become relatively more important as one moves away from the reso- 
nance. In either case, the LO, which ignores the range, contains energy-dependent errors 
of 0(/c^a2-R). In addition, in both cases, the truncation to a model space excludes physics 
of momenta beyond A, which introduces further energy- dependent errors of (9(/c^a2/A). In 
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other words, the truncation induces contributions to the effective range, the shape param- 
eter, and so on, governed by A rather than by 1/R. The role of contributions beyond LO 
is to correct for these two types of energy-dependent errors: NLO for k'^ errors and higher 
orders for higher powers of k. 

In the untrapped system, the power counting for a system with a large S'-wave scattering 
length 02 is such that corrections beyond LO have to be treated as perturbations j2[. It is 
one of our goals in this paper to show that treating higher orders in perturbation theory in 
the presence of the trap allows for a systematic improvement of the two-body energies. 

At NLO, we include corrections as first-order perturbations on top of the LO wavefunction 
As discussed above, in LO we fix Co so that one of the states has the "observed" 
energy. The NLO correction in Eq. ([3]) introduces a new parameter C2 that can be chosen 
so that a second energy level is fixed. However, the NLO term induces, in general, a non- 
vanishing correction to the energy used to fix Cq. One should, therefore, readjust Cq so 
that in NLO we reproduce the two observables (energy levels) at the same time. In order to 
keep track of this change, it is convenient to split Cq into an LO piece C'q \ which remains 
unchanged, and an NLO piece C'^\ as done in Eq. (fTO|) . and treat the latter in perturbation 
theory as well. 

Thus, for energies we have 



K 



(0)2 



C^o' + 4/iC 



K 



,.(0)2 



A^(0)2 



^(0)2 



(1) 



Ci^' + 4/iC 



(1) 



4d°) r 



'-^o 

3/253 



TT 



' I'Tnax 



n=0 



x + 



3%%^ r {rimax + 1) 



(24) 



where we used Eq. (1B15I) . 

The requirement that two energy levels have the correct positions in the spectrum fixes 
the amount of change from LO, thus providing two equations that determine the unknown 
coupling constants Cg^^ and €2^^ in each model space. In the case when the lowest level Eq 
is already fixed to a given (experimental or theoretical) value, Eq'^ = 0. However, in the 
case with finite non-negligible range, one can alternatively choose that in LO Cg^^ be fixed 
to a level of the two-body spectrum without a range, while in NLO that level is shifted to 
the correct position with range, thus requiring that Eq''^ 7^ 0. These two alternatives are the 
HO-basis equivalent to fixing Cg^^ in a free-particle basis to, respectively, a known binding 
energy (such as the deuteron binding energy) or the scattering length. In either case, if we 
take, say, the first excited level Ei{uj) to be reproduced at NLO in addition to the ground 



state, the two equations for the determination of Co^^(n^ 
written as 

Ei^\nmax,(^) = Ei{uj) - Ef^(n„a^,w), 
From Eqs. (l24l) and (l25l) . we can easily solve for CQ^\nmax 



. u) and C 



0,1. 



(1), 



n. 



max,uj) can be 



(25) 



, u) and C 



(1), 



4^ 

(-<(0)2 



(1) /,.(0)2 



.(0) 



(0) 



E^ 



(0) 



(26) 
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and 



(1) 



4G 



(0) 



r {ur, 



+ 



p(0) p(l) / (0)2 



(0) 



(1) 37r263 r 



+ 1) 



(1) /^(0)2 



(27) 



With these couphng constants fixed, Eq. fl24|) provides values for the other energy levels. 

The form of the NLO wavefunction is a bit more complicated; in App. Owe show that, 
up to higher-order terms. 



(0)/ 



(1)^ 



r,(0), 



E 



41/2) (^2/^2 



^0 J-(E(o)(n^,,) + E(i)(n^,,)) 



(28) 



where 



^(1) 



,(0)2 



TT' 



E 



W^(i50_^^ 



+ 2/i 



r^(0)2 



(29) 



Of course, just as in LO, the levels not used as input at NLO will have errors oc 1/A. 
As we show explicitly later, the magnitude of these errors is smaller than at LO, since more 
physics has been accounted for. If further precision is desired, we can continue the procedure 
to higher orders. 

In this paper we consider one order more, N^LO, so as to show the systematic trend of 
improvement — but we omit details that can be obtained straightforwardly, if painfully. The 
correction to the energy is obtained using perturbation theory up to the second order. 
In addition to the second-order correction from V^'^^'' (fTOj) . one has the first-order correction 
from y(2) (^^. 



E(2) 



(^(o)|y(2)|^(o)^ ^ 1 |(^(o)|y(i)|^(i)) + (^(i)|^(i)|^(o))} . 



(30) 



The potential V^'^^ affects the energy levels that were fixed already, so again it is convenient 

('2') 

to compensate for this by adding in Eq. (1T2|) perturbative shifts Cq and Cl with respect 

to the lower-order parameters. These two parameters, together with the four-derivative 

(2) 

parameter C4 , are determined so that three energy levels are fixed to the correct values. 



Taking the lowest three levels Ei{uj), 
determination of C, 



(2), 



0,1,2, to be fixed, the three equations for the 



and 



Ef\n^ax,^^) 



E^irir. 



0,1,2. 



(31) 



Obviously, higher corrections can be added in a similar fashion. 



Infinite-cutoff Hmit 



In the absence of a trap, the LO EFT is formally equivalent [2| in the infinite-cutoff limit 
to the pseudopotential [9|. As we show here, the situation is the same in the presence of the 
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HO potential, where the pseudopotential was solved in Ref . 10| . We consider explicitly the 
EFT to NLO, in order to derive also the first corrections to the pseudopotential in the trap. 

The wavefunction to NLO for a finite value of Umax is given in Eq. (128|) . By having 
n-max oo and using Eq. ( IB12p one obtains: 



+2y.«:W-|^5(r). (32) 

^0 



where U is the confluent hypergeometric function, and we have introduced the energy 

EW(oo,a;) + EW(oo,a;) 
e[oo) = ^ ^ (33) 

CO 

which is the limit of the energy (in units of u) of the two-body system in the trap. The 
second term in Eq. fl32l) was obtained by using the completeness of the HO basis, 

lim 0„,(O)0„(r) = lim ^ ^ (0|n/m)(n/m|f) = 5{r). (34) 

n=0 nl=0 m=—l 

The singularity of this term is mitigated by the pre-factor C2^V^o°^ which vanishes as 1/A 
for large cutoff. It stems from the enhancement of high virtual momenta due to the second 
derivatives in V^^\ As shown in App. |Al for an equivalent energy-dependent potential this 
term is absent. 

For small, non- vanishing values of r, i.e. < r ^ b, use of Eq. (1B14P gives 

^{0<r«b)^l{l- 2 ^): JJ, - + 0{-]}. (35) 

By identification of this wavefunction with the wavefunction in the untrapped case, Eq. (|2]), 
we obtain the relation between the energy of the trapped two-body system and the ERE 
parameters: 

^l--^£(oo) + ...|. (36) 




r (i - ^) 2a2 



In our framework, Eq. ( !36l) is to be interpreted in perturbation theory. At LO only the 
scattering length is taken into account, and the energies e^^\oo) are given by 



a relation first found in Ref. [10| using the pseudopotential [9[. The latter can be viewed 
as a renormalization of the delta- function interaction j2[. Indeed, the strength of the LO 
delta-function interaction has to decrease with the increase of the model-space size in order 
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for it to give sensible results. One can see from Eq. (IBlOp that for large n^ax 
in Eq. (123!) grows as ^Jrimax- Taking the limit of Eq. (123|) and using Eq. (|3^ 



(0) 



TT 



/iA 



7r 



2Aa2 



O 



( 1 



the first term 
, one finds 

(38) 



In the A6 — )■ oo limit this running is exactly the one found in free space [2| : the non-trivial 
rate of change with the cutoff A is controlled by 1 /a2 . This is not surprising since the large-A 
behavior should be independent of the long-range physics of the trap. 
At NLO, the effective range appears: Eq. (15^ becomes 



3 _ £(°)(oo)+£(^)(oo) 

2 



1 _ £(0)(oo)+£(l)(oo) 

2 



which can be solved to this order as 

2a2r2 



2a2 



^.(°)(oo) 



(39) 



r{0)| 



oo) 



62 



3 _ eW(oo) 

4 2 



e(0)(oo) >, ' 
4 2 



(40) 



in terms of the digamma function ?/;(°). In the presence of range, all levels change from LO 
to NLO. The range in fact controls the asymptotic behavior of the NLO coupling constants. 
In the large-rimai: limit we find from Eqs. ( 127|1 and (126|) that 



7r^r2 
'T2)7 



1 + C 



/ 1 



1 



Vr2A' a2A 



and 



/iy2 
Stt 



r2A 



(41) 



(42) 



Again, this is the same running as in free space |2[, as it should be. Note that it qualitatively 
changes for r2 = 0. In this case, e'^^^[po) = 0, which means, given our choice of levels to 



fix the coupling constants at NLO, that Eq'^ = and i?}^"* = 0{Ef/A). As a consequence 



(1) 



the ratio of energies in Eq. fl2Bl) goes to as A — )■ oo, and /xACg and yuA^C2 approach 
constants. 

At N2L0 the power counting of Ref. ^ indicates that no new ERE term should be 
included. This is a reflection of the fact that a fine-tuning in 02 does not in general lead 
to an enhancement in the shape parameter P2: in order to go from r2/c^/2 to P2k^/^ two 
orders are needed. Therefore at N2L0 Eq. (!36l) is unmodified. It is only at N^LO that we 
need to account for a non-vanishing S'-wave shape parameter, also with an S'-wave potential 
of the form of Eq. (IT^ (but with different parameters). For fermions, N^LO also contains 
a P-wave interaction to account for the P-wave scattering volume. 

It is clear that the procedure can be continued to higher orders, and we expect for / = 
levels 



2 



hk 



cot 5o(^), 



(43) 
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where 

kh = ^2e{oo) (44) 
and 5n(/c) is given by the ERE ([T]). This extension to subleading orders agrees with Refs. 



llHl3|. Waves with / > 1 can also be examined with the same method we developed here, 



but we leave details for future work. 

The importance of Eq. ( 136|) lies in the link between the energies inside the HO well, £^(oo), 
and the scattering parameters, 02, ^"2, etc. It is the analog of Liischer's formula [l6[, which 
links the levels inside a cubic box and the same scattering parameters. As such, Eq. ( 136|1 
provides the energy levels necessary to fix the coupling constants in finite model spaces when 
the ERE parameters are known. In the nuclear case, for example, the scattering parameters 
have been determined from data, so Eq. f l36p can be used to fix theparameters of the 
pionless EFT without relying on a fit to nuclear levels, as done in Ref . [7| . An extension of 
Ref. 0] to this case is in progress. 

In the atomic case, the lowest energy Eq{uj) of two trapped particles has been measured 
[ij. One can use Eq{uj) directly as input in Eq. f l23|) . Alternatively, this energy was found [l[ 
to be in good agreement with the lowest state of the theoretical spectrum obtained under a 
pseudopotential assumption [lo[, in which the eigenvalues are determined by the scattering 
length a2- In Ref. [sl we used this theoretical LO energy as input; we confirmed that the two- 
fermion spectrum of the underlying pseudopotential is reached asymptotically as rimax 00, 
and we calculated the properties of three- and four-fermion systems. In subleading orders, 
we can now include effective range and higher effects, which might account for the small 
discrepancies [l| between the theory and experiment. 

We should stress that, because the bare coupling constants are not observables, their sizes 
have no direct physical meaning. What matters is the total contribution of a given order to 
an observable. For example, for A> l/r2, Cq°^ ^Cq^'', and yet, the NLO energy shift ( 140|1 
is small as long as r2 is sufficiently small. We show explicit results for energies in the next 
section. 



III. RESULTS 

In this section, we illustrate the approach described above in a few cases of interest 
for atomic and nuclear systems. We initially consider situations where the range of the 
interaction can be neglected; we first look at the unitarity limit, h/a2 = 0, and then finite 
b/a2- Finally, we consider the case where the range of the interaction is non- negligible albeit 
still small with respect to the scattering length. In all cases, we use Eq. fl5B|) to provide the 
asymptotic levels in the trap. 



A. Unitarity limit 

We start by considering a system of two trapped particles at unitarity, characterized by 
b/a2 = and r2/b = 0. In the untrapped case, this situation can be realized by consid- 
ering an attractive potential given by, for instance, a square well with a fine-tuned depth. 
Asymptotically, the wavefunction at zero energy behaves as 1/r and the scattering length 
a2 is then infinite. 

In the presence of the harmonic trap, the only scale is set by b, or, equivalently, u, so 
the solutions of Eq. (I36l) have to be constants; indeed, they are given by the poles of the 
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Gamma functions in the denominators, i.e., 

en(oo) = ^ + 2n, (45) 

where n > is an integer. At each order we use a finite number of these energies to determine 
the interaction parameters in each model space: LO, Sq and ei at NLO, and Eq, £i, 

and £2 at N^LO. 

Using the ground-state energy to determine Cq"'*, Eq. f l23|) simphfies to 

^(0). X TT^fe r {n^a^ + 1) 

Q (nmax,CU) = -y, (46) 



which at large cutoff becomes 



2/ir 



i + o 



(47) 



Analogous expressions can be derived for the other coupling constants. The resulting running 
of the coupling constants is shown in Figs. [H |21 and El One can see that, in agreement with 
Sec. Ill C[ a coupling constant Cj behaves asymptotically as 1/A^*+^. 

With the interaction parameters so determined, we calculate the remaining energies, 
which depend on rimax- Results for some excited states, from the second (n = 2) to the 
fourth (n = 4), are shown in Fig. IHfor different values of the cutoff rimax in the dimensionless 
combination Ab. At LO and NLO, all states change with Ab; at N^LO, the second excited- 
state energy is used as input. As we showed in Sec. Ill CI these calculated energies converge 
as Umax increases to the values given in Eq. (145|) . The plot explicitly shows, in addition, 
that the convergence with respect to rimax to the exact value is increased as higher-order 
corrections are added in perturbation theory. 

In order to further study the effects of a finite rimax, we consider Eq. (15^]) from a different 
angle. Hitherto, we have used it simply to fix observables (energy levels). Now we employ 
it also to extract phase shifts. In a given model space, characterized by a given rimax, we 
calculate the energy levels en{rimax)- The lowest levels used in the fitting of coupling con- 
stants are, of course, exact, while all others deviate from the exact values. The momentum 
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10 15 20 

Ab 



FIG. 2: The coupling constants fiA^C^^ and fiA^C^^ at unitarity as a function of the cutoff in the 
dimensionless combination A6. 




FIG. 3: The couphng constant /lA^C^ 
combination A6. 



(2) 



at unitarity as function of the cutoff in the dimensionless 



kn associated with an energy level Er 



[Ur. 



IS 



kr,h 



(48) 



so that the phase shifts for discrete values of /c, 5o(^n.), can be determined by simply inverting 
Eq. (USD: 

1 h ^Kii \ ^ r(3/4-£„(n^a^)/2) 

knh cot 5o[kn) = -2—— — — -. 49) 

r(l/4 - en[nmax)l2) 

This allows us to study the effect of the truncation in a HO basis in terms of the more 
familiar ERE parameters: any deviation from the value fcfecot (5o(^) = is due to truncation 
errors and we can calculate the induced range, shape parameter, and so on. In Fig. O we 
plot kb cot Solk) in a model space with Umax = 15 as a function of k'^b'^. (It is more natural 
to express k in units of l/a2, but at unitarity b provides the only length unit.) At LO, 
kb cot Soi^k) starts off linear in k^b^: a linear fit (indicated by the dashed line) shows an 
induced effective range of about 0.176 in this particular model space. At larger values of 
k'^b'^ deviation from the linear behavior is seen, indicating the presence of further induced 
ERE parameters. NLO and N^LO corrections reduce the size of the ERE parameters, so 
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Ab Ab Ab 

FIG. 4: Running with the dimensionless cutoff Ab of the second (£2)5 third {£3), and fourth (£4) 
excited-state energies, in units of w, for two particles in a trap with the interaction characterized 
by 6/02 = and r2/b = 0. Results at LO (green diamonds), NLO (red squares), and N^LO (black 
circles) are compared with the exact values (blue solid lines) given by Eq. (f36l) . The coupling 
constants are fixed so that the ground-state, additionally the first excited-state, and additionally 
the second excited-state energies are reproduced in LO, NLO, and N^LO, respectively. (Note that, 
indeed, at N^LO the second excited state in the leftmost panel is constant on top of the exact value 
for all values of Ab.) 

that the results for the phase shifts improve order by order, getting closer and closer to the 
horizontal axis. Since here A%^ = 63, at ~ 60, the errors are dominated by the higher 
orders, and therefore the lower orders do little to improve on the previous orders. However, 
at low momentum, the results systematically improve, as expected. 

Our perturbative treatment thus provides a way to systematically reduce the effect of 
truncation to a reduced model space, demanded in the calculation of larger systems. Em- 
boldened by this success, one might be tempted to treat the subleading potentials exactly in 
the Schrodinger equation. Apart from the renormalization problems pointed out in Ref . 0] , 
we find no obvious numerical improvement when this is done here. In Fig. also show for 
comparison results obtained when the NLO correction to the potential is not considered as 
a perturbation but fully diagonalized together with the LO potential. (The same two lowest 
levels are used to fix the two coupling constants of the LO+NLO potential.) One can see 
that by considering the LO+NLO potential this way, results are further away from the exact 
curve than by treating the NLO potential as a perturbation. This result is not particular 
to this example; more generally, truncation errors get worse when subleading corrections 
are treated improperly. The reason is that doing so includes only part of the higher-order 
corrections; it neglects the rest, needed to ensure systematic improvement. 

B. Finite scattering length and zero range 

In the previous subsection we saw that the rate of approach to the asymptotic values of 
two-body energies improves systematically as the order increases at unitarity. We now show 
that there are no qualitative changes when we consider the case of a pseudopotential away 
from unitarity [lo|, where the scattering length 02 is finite, with r2 still vanishing. 
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FIG. 5: S'-wave scattering phase shifts kb cot 5o{k) as function of the dimensionless squared mo- 
mentum k'^b'^ in the finite model space characterized by Umax = 15, at unitarity. The points at 
LO (green circles), NLO (red squares), and N^LO (black diamonds) are obtained from calculated 
energies via Eq. (j49|) . The dashed line corresponds to a linear fit of the LO curve at small kb 
values. For comparison, we also show results with the NLO potential fully diagonalized (black 
stars). 

The corrections to the potential are taken into account as in the unitarity case, i.e., the 
LO correction is iterated to all orders whereas higher corrections are treated as perturbations. 
The parameters at each order are adjusted so that the lowest levels satisfy Eq. exactly. 
The A dependence of coupling constants is similar to unitarity, except for a markedly slower 
convergence for Cq. This is particularly obvious for Cq^\ when we compare the more general 
Eq. f p8|) . which applies here, with its unitarity version f H7|) . 

As for energies, let us first consider a2 < 0. This situation corresponds to the case when 
the depth of the potential between the two particles is decreased starting from the fine-tuned 
value at unitarity. As an example, we show in Fig. |6] results for b/a2 = —1 for the same 
energy levels previously displayed at unitarity in Fig. |H The exact values are slightly higher 
than at unitarity. As before, convergence of the energies to the exact values is improved as 
more corrections to the potential are added, and the difference between the truncated-space 
energy and the exact result is mitigated as more corrections are included. 

For a weak enough potential, \a2\ becomes small. The energies are then close to the HO 
energies, 

3 

en = - + 2n + fe„, (50) 

where fe„ is a small correction and n = 0, 1, . . .. By using Eq. fl36l) and expanding the 
Gamma function around its poles one finds 



where 
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FIG. 6: Same as in Fig. [H but for 6/02 = — 1- 




FIG. 7: Same as in Fig. HJ but for h/a2 = 1. 



We have verified numerically that in this case convergence to the exact value can even be 
sped up by considering all interactions as perturbations, as in the "natural" continuum case 
discussed in Ref. [I]. 

We now turn to the case where 02 > 0. As the interaction between particles becomes 
stronger (starting from the case at unitarity) the scattering length 02 decreases. The ground 
state can have negative energy. The plot of the energy for excited states (starting from the 
second) for h/a2 = 1 is shown in Fig. [71 The exact values are now slightly lower than at 
unitarity. Again, there is no qualitative change in the pattern of convergence with respect 
to unitarity or negative 02- 

For a strong enough interaction, the absolute value of the ground-state energy, |eo(oo)|, 
is large and in the trapped system we have 



r(3/4-eo(oo)/2) 
r(l/4-£o(oo)/2) 

Using Eq. (15^ . we find 

£0(00) = 



£o{oo) 



O e 



00, 



162 
'2^ 



(53) 



(54) 
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Thus, when 6/02 — > 00, the ground-state energy Eo^rimax, ^) converges to £"0(00) = — l/2/ia|, 
the untrapped bound-state energy for the case 02 ^ R. The excited levels are again close 
to the HO values 

3 _ . 1 



2+2(^-1; 



-- + 2n + 5er 



(55) 



where n = 1,2, 
corrections 



Following the same procedure as for large, negative b/a2, we obtain the 

4 02 



>i 



r(|-n)p„^ 



1 + 



(56) 



C. Interaction with finite range 

A finite interaction range R usually generates higher ERE parameters of the same mag- 
nitude, \r2\ ~ R, IP2I ~ R^, etc., even when there is fine-tuning that leads to \a2\ ^ R. As 
a last example, we account for a finite effective range r2. 

Regardless of the quantity chosen as input in LO, the existence of range introduces errors 
that are energy dependent and can only be accounted for in subleading orders. As before, 
we use as LO input the ground-state energy Eq given by Eq. f l371) . This fixes the running 
of Ci°^ to the same values as in the previous subsection. However, at NLO, we obtain Cg^'' 
and 6*2^^ from the first two states of Eq. fl36|l with a2 and r2 non-vanishing. NLO is, thus, 
different from the previous subsection: it accounts not only for errors 0{a2k^/A) due to the 
explicit truncation to the model space but also for implicit ones, 0{a2Rk'^), in the potential. 

As discussed in Sec. Ill C[ the N^LO corrections are slightly more subtle 0]. The intro- 
duction of range leaves an error that can be as big as Ola^r'^k'^) = 0{a\R?k'^). Errors from 
the explicit truncation of the model space are now 0{a\k^ / K^) or 0{a\r2k'^ / and, as in 
general, are smaller than errors from the truncation of the expansion once A > These 
types of errors are one order in kR or k/A from NLO, which requires V^"^^ for control. In 
contrast, errors from the shape parameter are only 0{a2P^k'^) = 0{a2R^k^), two orders in 
kR down from NLO. Thus, at N^LO we determine Cq^\ C2^\ and from the lowest three 
levels of Eq. fl36|) . still with non- vanishing 02 and r2 and neglecting all higher-order ERE 
parameters. 

As an illustration, we take r2/b = 0.1. The coupling constants for finite scattering length 

a2/b = 1 are plotted in Figs. [SI [HI and[TUJ The running of Cq'^ is identical to the one in the 
previous subsection, and as advertised is gentler than the one displayed in Fig. [TJ As seen 
in Eqs. f HT]) and f H2|) . the range changes the running of Cq^^ and 02^^ dramatically with 
respect to the zero-range case, say Figs. [1] and |2l The new runnings approach the limits 
given in Eqs. (jH]) and (gl]). At Ah ~ 38 we find, for example, SttC^^^ / j2r2C^°^ ~ 0.7. The 
slow convergence is a consequence that at this cutoff r2A is only ^ 3.8, and the 0{l/r2A) 
corrections are still not so small. For the other coupling constants, Cq , C2 , and C4 , 
the running is also very different from the zero-range case. In agreement with our error 
estimates, these parameters are nearly flat when normalized with l/r"^. 

Qualitatively, the only change in energies with respect to previous subsections is the finite 
jump from LO to NLO due to the effective range. Energies for a few excited states in the case 
of an infinite scattering length, b/a2 = 0, are plotted in Fig. [TTl and for a finite scattering 
length, b/a2 = 1, in Fig. [T2j As evident, LO misses the correct asymptotic behavior by 
a little bit because it lacks the effective range. The NLO results converge, as rimax — ^ 00, 
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FIG. 8: The coupling constants ;uACg'^\ fiC^^^ /r2, and /xCQ^^/r^A for 6/02 = 1 and r2/b = 0.1 as 
function of the cutoff in the dimensionless combination Ab. 
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FIG. 9: The couphng constant ^Is?C2^ /r2 and iiAC^^ jr^ for 6/02 
of the cutoff in the dimensionless combination Ab. 



1 and r2/6 = 0.1 as a function 



to the values given by Eq. (15^ with the range, as they should. However, the inclusion of 
N^LO corrections speeds up convergence considerably: for a fixed value of the cutoff Umax 
the results at N^LO are much closer to the exact value. It is straightforward, for example, 
to generalize Eqs. f l5^ and fl56|) to non-zero range. 



IV. CONCLUSIONS AND OUTLOOK 

We have extended the work initiated in Ref. jsf on the application of effective field theory 
to the trapped two-particle system to include subleading orders. In EFT, the short-range 
interparticle potential is replaced by a delta function and its derivatives, as we presented 
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FIG. 11: Same as in Fig. [J but for r2/b = 0.1. 



in the main text, or equivalently by an energy- dependent delta function, as discussed in 
App. [Al The singular nature of the interaction requires regularization and renormalization, 
which is naturally accomplished by the use of a finite model space, as required in actual 
calculations. Once the parameters of the interaction are fixed in each model space using a 
finite number of levels, other levels can be calculated. 

The leading-order interaction is solved for exactly, whereas higher contributions are 
treated in perturbation theory. We have considered the corrections to the potential up 
to next-to-next-to-leading-order. We have shown explicitly that this method can systemati- 
cally account for the physics of the effective-range expansion, while treating the subleading 
potential exactly gives worse results. 

In the limit of a large model space, the leading-order theory reproduces the pseudopo- 



tential result of Ref. |10| . while in subleading orders our approach gives range and higher 
ERE corrections lll-ll3| to the pseudopotential result. Therefore we have provided an alter- 
native derivation of these results. In turn, these asymptotic results allow us to incorporate 
the physics of scattering into the trapped system: by fitting the coupling constants to the 
asymptotic values of some levels, we can calculate other energies without necessarily resort- 
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FIG. 12: Same as in Fig. HJ but for h/a2 = 1 and r2/b = 0.1. 



ing to fitting measured bound-state energies. This is important in the nuclear case where 
the scattering parameters are known. 

We have studied numerically the convergence to large model spaces in some detail. We 
have presented results at unitarity and for finite values of the scattering length 02 as well 
as for finite values of the effective range r2- In all cases, we have observed convergence to 
the asymptotic values. Moreover, for a fixed value of the cutoff, the difference between the 
exact value and the value obtained in a finite model space is mitigated as more corrections 
to the potential are taken into account. Truncation errors can thus be reduced not only by 
increasing the size of the model space but also by increasing the order of the calculation. 

Thus, we have a method to calculate the energies of few-body systems that is independent 
of the details of the short-range potential and can be improved systematically. Although 
in the two-body system our method is simply an implementation of the effective-range 
expansion, it is well-suited for extension to larger systems where numerical calculations 
are required. In these cases, limited computational power restricts the size of accessible 
model spaces. We can fight this limitation by increasing the order of the calculation, thus 
accelerating convergence. Results for systems with more particles will be presented in a 
future publication |14| . 
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Appendix A: Potential 



Analogously to the multipole expansion in electrodynamics, the short-range two-body 
potential can be expanded in a power series in momenta j2|, 17| . 



+C'^p' ■ p + C',p' ■ p(p'^ + p^) + Cl{p' ■ + . . . (Al) 

where p (p ') is the initial (final) relative momentum, the Cf are constants, and ". . ." denote 
terms with six or more powers of momenta. Here the Cq, C2, and C4 terms contribute to the 
5* wave, C2 and C4 to the P wave, and C'l to the D wave. The contribution from to the 
on-shell two-body system vanishes, and thus it is only relevant in larger systems, where it 
cannot be separated from few-body forces. If we were to include this term and set a similar 
linear system to Eq. fl3T]) (based on the fit of the first four levels inside the HO), this system 
would have no solution: its zero determinant would indicate that the corrections C4 and C4 
are not independent. 

In the situation of interest here, where the two-body 5- wave scattering length 02 is large, 
the S'-wave constants are enhanced over the others by powers of a2, and as a result the 
Co term is LO, the C2 term is NLO, the C4 term is N^LO, while the others appear only at 
higher orders |2|. In this paper we limit ourselves to N^LO, where only the S wave is present. 
Taking the Fourier-transform with respect to both p and p\ the potential is found to be 
given in coordinate space by Eq. ([3]). It is non-local in the sense of involving derivatives. 

A completely equivalent formulation of the EFT in the two-body system is achieved with 
an energy- dependent potential [2|. It can be implemented through an auxiliary "dimeron" 
field [18|. In LO, the dimeron is characterized by a mass A and a coupling constant (7 to a 
two-particle S'-wave state. In the case of interest here, |a2| ^ -R, the kinetic energy of the 
dimeron is an NLO effect j2|. A subtlety is that the bare dimeron can be a ghost: the sign a 
of the kinetic term can be positive or negative, depending on the sign of the effective range. 
Denoting by E is the total energy in the center-of-mass frame, the two-body potential is 
simply 



V{E) = + (A2) 



in momentum space, and 



V{f',r) = fl + + . . . ) SinSir) (A3) 



A V A 

in coordinate space. It is local, but energy dependent. 



This potential was considered in Ref. [13[. An alternative treatment follows the same 
steps as that of the momentum-dependent potential in the main text, with the substi- 
tutions 

,(0)2 

AW' 



- (A4) 



A(i) V AW, 
Ci'^ (A6) 
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In particular, the LO is identical to that presented in the main text, with Eq. f lA4p . Note 
that since only the ratio ^((''^^/A^^^ enters, the separation between g^^^ and A*^*^^ is arbitrary. 
At NLO, one finds equations that are somewhat simpler than Eq. (12^ and the ones that 
follow. With Eqs. (lASp and (IA6p one can find the renormalization of g^^^ and A^^^ separately. 
The appearance of energy instead of momentum leads to a softening of the UV behavior. 
This results in a wavefunction at NLO of exactly the form fl28l) and (1291) . but with Eq. 
(lA6p . In the infinite-cutoff limit the most singular term in Eq. fl5^ is absent, and Eq. fl35p 
follows directly. Any observable in the two-body system is identical for this potential and 
the momentum-dependent potential discussed in the text. 



Appendix B: HO notation and definitions 

The HO basis functions with the length parameter 

1 



are the solutions of the three-dimensional Schrodinger equation 



(Bl) 



1 

2 



with energy (in units of u) 



They are given by 



^ ^ + 52 J ^nlm{r) = enl(f)nlm.{r) 



Enl o I 7 I ^ 

eni = = 2n + / + -. 

uj 2 



(B2) 



(B3) 



{f\nlm) = 4>nim{r) = Rni{r)Yim{r), 



(B4) 

where Yim{r) are the usual spherical harmonics, and the radial parts Rniij) can be shown 
to have the form 



1/2 



Rni{r) 



^63r(/ + 3/2)^ 



1/2 /r\' 



4'+V2) (0)] ' M exp (-rV262) 4'+V2) (^2/^2^ ^^35) 



The Z/^")(x)'s are the generalized Laguerre polynomials, which can be written as 19 

Tin + a + 1 



Li") I 



X) 



n\ r(a + 11 



-M{-n,a + l,x) 



(B6) 



in terms of the confluent hypergeometric function M. 

For contact interactions, it is useful to know the value of the radial wavefunction at the 
origin. Because of the r' factor, only / = contributes, as expected. The S wavefunction of 
energy En = {2n + 3/2)a; is, omitting the I = m = labels. 



(f |n) = Mr) = 7r-'/'r3/2 [4^/2) (0) 
In particular. 



-1/2 



exp 



^2/252\ ^1/2) / 2/^2\ 



MO) = TT-'/'b-'/' fL(V2) (0) 



1/2 



(B7) 



fB8) 
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where 



4°'(o)- ^„':v:.°v,' =n(i+^ 



n!r(a+l) 

Use of Stirling's formula yields 



i + -)4"i(o). 

n' 



r(z + a) 
for large which leads to 



2" <^ 1 



aioL — 1) 
2z 



77, 



r(a + 1) 



2?7, 



(B9) 



(BIO) 



(Bll) 



for large n. 

The generalized Laguerre polynomials satisfy [lo| 

E^^^ = r(«)f/K 3/2,0;) 



n=0 



n + a 



in terms of the confluent hypergeometric function pi 

(a- 1/2, 1/2, x) ViM (a, 3/2, x) 



?7(a,3/2,a;) 



For small x, 



r(a)f/(a,3/2,x) 



r(a)r(l/2) r(a- l/2)r(3/2) 

r(a) 



(B12) 



(B13) 



1 - 2- 



X + 0{x) 



r(a- 1/2) 

Also useful are sums involving the generalized Laguerre polynomials at the origin 

n=0 



(B14) 



E 



" ' ' 3v^ r(m + l) 
41/2) (0) 4 r(m + 3/2) 



n + 1/2 r(m + l) 



(B15) 
(B16) 



2v/^ 



r(a) 



;^o(n + l/2)(n + a) 

r(m + 3/2) 3F2 (1, m + 3/2, a + m + 1; m + 2, a + m + 2; 1) 



r(a+ 1/2) 



and 



r(m + 2) 
L(V2)(0) 



-2^ 



\ r(a) 



7r(a + m + 1) 



(BIT) 



S(^ + l/2)(n + a)2 ^^"\r(a + l/2) 

r(m + 3/2) 4F3 (1, m + 3/2, a + m + 1, a + m + 1; m + 2, a + m + 2, a + m + 2; 1) 



r(m + 2) 



7r(a + m + 1)2 



(B18) 

where and 4F3 are generalized hypergeometric functions and "^^^^ = F'/F is the digamma 



function 19 
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Appendix C: Wavefunction at NLO 



The correction at NLO is taken into account as a perturbation to the LO potential. We 
write the corresponding quantum state of the two-fermion system as = \ip^^^) + 14^^^^) 
where \^^^^) is the solution at LO and IV'*-^-') the correction at NLO. Both of these vectors 
are expanded in a HO basis, as in Eq. f lT^ . 

Corrections at NLO are obtained by solving the Schrodinger equation at first order, Eq. 
( ITT]) . Let us introduce the projection i?„ of the left and right sides of Eq. (fTT]) . 

Rn = - = - (CI) 

and define k'-^-* as 

1) = C^^V^^HO)- From the Ihs of Eq. ([CT]) we can write the coefficient 
of the wave function at NLO as 



Normalizing the wavefunction at LO to unity, 

0^(0) 



K 



(0)-2 



E 



{En - Ewy 



(C2) 



(C3) 



Keeping this normalization at NLO, l^^'^^) and l^/''-^^) are orthogonal; one then obtains, from 
Eqs. ([HD for and ^ for c^, 



(0)„(1) 
n n 



,(1) 



K 



-(0) 



E 



^>n(0)i?„ 



0. 



^ {En - Ei^)y 

On the other hand, the rhs of Eq. ( ICll) provides an expression for Using Eq. ([18 



(C4) 



[VV^°Hr)],=o = E„0„(O) = 2/i 



Inserting this expression into the rhs of Eq. fICip and eliminating tjj^^\0), 



Rn = 0n(O) 



(0) 



(0) 



y - C« + 2/.C« (En + - E ^UO) 



(C5) 



(C6) 



From Eqs. (lC4p and (]C6p . we arrive at an expression for k*^^^ using Eq. (120|) : 



(1) 



,(0)3 



0^(0) 



(1) 



E 



;(o) 



+ 2/iC; 



(1) 



,(0) 



(0) 



On n 



(C7) 



We can now finally obtain the coefficient c^^-* by plugging Eqs. f lC6p and f lC7p into Eq. (\C2\i : 



.(1) 



.(0) 



^„ - Eio) 



(0) 



(i?(o) - E„ 



2/iC 



(1) / n 



(0) 



-2/i«:W-^0„(O). 

^0 



(C8) 
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The total wavefunction at NLO is then given by 

n 



(C9) 



n \ 



a 



(1) 



where A^^^ is a first-order term defined as 



r) + 2/iKW-^E0„(O)0.(r), (CIO) 



^(1) = ^(0)2 



.(0) 



2/i 



(0)2 



(Cll) 



At first order in A^^^^ and E'^^), the previous expression for the wavefunction is equal to 



0n(O)0„(r) 



(1) 



^(0) + ^(1) ' MO) 



+ 2^K^'^-^Y.M0)Ur). (C12) 



Upon insertion of the expression for the radial wavefunction of the HO basis we obtain Eqs. 
(1^ and dsn]). 
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